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Density functional theory (DFT) embedding provides a formally exact framework for interfacing correlated 
wave-function theory (WFT) methods with lower-level descriptions of electronic structure. Here, we report 
techniques to improve the accuracy and stability of WFT-in-DFT embedding calculations. In particular, 
we develop spin-dependent embedding potentials in both restricted and unrestricted orbital formulations 
to enable WFT-in-DFT embedding for open-shell systems, and we develop an orbital-occupation-freezing 
technique to improve the convergence of optimized effective potential (OEP) calculations that arise in the 
evaluation of the embedding potential. The new techniques are demonstrated in applications to the van-der- 
Waals-bound ethylene-propylene dimer and to the hexaaquairon(II) transition-metal cation. Calculation of 
the dissociation curve for the ethylene-propylene dimer reveals that WFT-in-DFT embedding reproduces full 
CCSD(T) energies to within 0.1 kcal/mol at all distances, eliminating errors in the dispersion interactions due 
to conventional exchange-correlation (XC) functionals while simultaneously avoiding errors due to subsystem 
partitioning across covalent bonds. Application of WFT-in-DFT embedding to the calculation of the low- 
spin/high-spin splitting energy in the hexaaquairon(II) cation reveals that the majority of the dependence 
on the DFT XC functional can be eliminated by treating only the single transition-metal atom at the WFT 
level; furthermore, these calculations demonstrate the substantial effects of open-shell contributions to the 
embedding potential, and they suggest that restricted open-shell WFT-in-DFT embedding provides better 
accuracy than unrestricted open-shell WFT-in-DFT embedding due to the removal of spin contamination. 



I. INTRODUCTION 

The demand for accurate and efficient descriptions 
of complex molecular systems requires development of 
quantum embedding methods for electronic structure in 
which a small subsystem is treated with a high level 
of theory while the remainder of the system is treated 
at a more affordable level. Widely used examples of 
quantum embedding include QM/MMjD® ONIOMpsl 
and fragment molecular orbital (FMO) approaches J2HE0 
which have led to significant advances in the simulation 
of condensed-phase and biomolecular systems. However, 
such methods generally rely on empirical models for the 
subsystem interactions, including link-atom approxima- 
tions for embedding across covalent bondiP^HHIancl point- 
charge electrostatic descriptions of the environment,^! 
that are difficult to systema tically im prove and that can 
fail in practical applications ! 3 * 6 * 16 -*^ 

Density functional theory (DFT) offers an appealing 
framework for addressing this challenge EHEH DFT em- 
bedding provides a formulation of electronic structure 
theory in which subsystem interactions depend only on 
their electronic densities, including non-additive contri- 
butions due to the electrostatic, exchange-correlation 
(XC), and kinetic energy terms. In the WFT-in-DFT 
embedding approach, the DFT embedding potential is 
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included as an external potential for WFT calculations, 
providing a WFT-level description for one (or more) sub- 
system while the remaining subsystems and their inter- 
actions are seamlessly treated at the DFT level of theory. 

Several groups, including this one, have recently 
demonstrated that non-additive kinetic energy contri- 
butions to the em bedding potential can be exactly 
computecP^HESM] t ne use f optimized effective 
potential (OEP) methodsp2lEU In this paper, we intro- 
duce a simple technique to improve the robustness of 
OEP calculations in systems that exhibit small HOMO- 
LUMO gaps, such as transition metal complexes. In ad- 
dition, we derive spin-dependent embedding potentials 
to enable the accurate description of open-shell systems 
in the WFT-in-DFT embedding framework. Numeri- 
cal applications to the van-der-Waals-bound ethylene- 
propylene dimer and to the hexaaquairon(II) transition- 
metal cation illustrate the applicability of these new tech- 
niques and demonstrate the accuracy of the WFT-in- 
DFT approach in systems for which conventional density 
functional theory methods exhibit substantial errors. 



II. THEORY 

Like Kohn-Sham (KS)-DFT, DFT embedding provides 
a formally exact framework for the ground-state elec- 
tronic structure problem. Here, we review DFT-in-DFT 
embedding and its basis for WFT-in-DFT calculations. 
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A. DFT-in-DFT embedding 

We begin by considering a closed-shell system in which 
the total electronic density pab consists of two subsys- 
tems, pab = PA + Pb- The corresponding one-electron 
orbitals for pa and pe obey the Kohn-Sham Equations 
with Constrained Electron Density (KSCED)j22 
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where i = I,..., N A , j = 1, . . . , iV B , and N A and N B 
are the number of electrons in the respective subsystems. 
v c e is the effective potential for the coupled one-electron 
equations, 



"cfi [pa, Pab; r] = v cS [p A ; r] + u e mb(r), 



(3) 



where v^[pa', r] is the standard KS potential for subsys- 
tem A, and 

«emb(r) = ^icM + v -Apb] r] + v xc [p AB ; r] - 

v xc [p A ; r] + u nad [p A , pab ; r] . (4) 

Here, i> B e (r) is the nuclear-electron Coulomb potential 
from the nuclei contained in subsystem B, vj is the 
Hartree potential, and v xc is the XC potential. In ad- 
dition to these familiar terms from conventional KS- 
DFT calculations, DFT embedding introduces the non- 
additive kinetic potential (NAKP) which properly en- 
forces Pauli exclusion between the subsystem densities. 
It is obtained from 
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where T s llad [p A , p B ] = T s [p AB ]-T s [p A ]-T s [p B }, and T s [p] 
is the non-interacting kinetic energy functional. The to- 
tal energy functional for the full system is then 



E[ P ab] = T s [p A } + T s [p B ] + T s nad [p A , Pb] 
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where the last three terms on the right-hand side (RHS) 
are the nuclear-electron Coulomb energy, the Hartree en- 
ergy, and the XC energy computed over the total density. 
Enforcing w em b(r) to be identical for all subsystems (see 
Sec. Ill B) leads to a unique partitioning in the DFT em- 
bedding formulation, such that the specification of the 
nuclei and the integer number of electrons in subsy stem 
A and B fully specify the density partitioning. 26 35 Eqs. 
1-|6] are easily generalized to the description of multiple 
embedded subsystems. 

We have previously demonstrated that by using OEP 
methods to calculate the NAKP, DFT-in-DFT embed- 
ding can accurately describe both weakly and strongly 
interacting subsyst ems, including subsystems connected 
by covalent bonds EZES1 anc j we have shown that this 



method is computationally feasible for large systems by 
way of localized approximations to the NAKPP^I More 
recently, we have introduced a projection approach that 
completely avoids the NAKP calculation in exact DFT 
embedding} 2 ^ which appears worthy of further investiga- 
tion. The OEP-based approach employed here is appeal- 
ing because it provides a local embedding potential that 
is a functional of only the subsystem electronic densities. 

In practice, the KSCED equations (Eq. [T] and Eq. [2]) 
are solved by simply modifying the core Hamiltonian in 
the self-consistent field (SCF) calculation to include the 
additional embedding terms. The embedding potential 
(Eq. [1]) can be written in the atomic orbital (AO) basis 
as 
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where the various terms on the RHS of this expression 
correspond to those in Eq. [31 The subsystem and total 
AO density matrices in Eq. |7j satisfy 7a + 7b = 7ab • 
It follows that the Fock matrix for subsystem A can be 
expressed as 
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where 
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and h A is the core Hamiltonian for subsystem A (the 
kinetic energy plus external potential due to the nuclei in 
subsystem A). The Fock matrix for subsystem B, f B m A , 
is analogously defined. 



B. WFT-in-DFT embedding 

The embedding potential in Eq. [4] describes the sub- 
system interactions in terms of their corresponding elec- 
tronic densities. However, the subsystem densities can 
be computed with any level of theory, thus allowing for 
the description of one subsystem at the (single- or multi- 
reference) WFT level, while the remaining environment 
is treated at the DFT level! 23 | 29 | 34 | ^ nH5l closed-shell 
WFT-in-DFT embedding simply involves performing a 
WFT calculation on a given subsystem using the modi- 
fied core Hamiltonian, h A 111 B in Eq.|9] that contains the 
embedding terms due to the environment of the other 
subsystem. The WFT-in-DFT energy is then obtained 
by modifying the DFT energy with respect to subsystem 
contributions at the WFT levelj^ 
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This expression is easy to evaluate since the terms in 
the parentheses are just the DFT and WFT energies of 
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subsystem A performed using the modified core Hamil- 
tonian, h AinB . Just as DFT-in-DFT embedding is an 
exact theory for the case of an exact DFT XC functional, 
WFT-in-DFT embedding is an exact theory for the case 
of an exact DFT XC functional and a full configuration 
interaction (FCI) WFT description!^! 

III. METHODS OF IMPLEMENTATION 

Here, we describe techniques to improve the accuracy 
and convergence of both DFT-in-DFT and WFT-in-DFT 
calculations. First, a description of open-shell DFT em- 
bedding is developed to incorporate the effects of spin- 
dependence in the embedding potential. Then, imple- 
mentation of the OEP calculation is discussed, and an 
orbital occupation constraint is introduced to enable ro- 
bust DFT-in-DFT and WFT-in-DFT embedding calcu- 
lations for systems with low-lying virtual orbitals, such 
as transition metal complexes. 

A. Embedding for open-shell systems 

For an open-shell embedded subsystem, the a and /? 
electrons generally experience different embedding poten- 
tials due to differing non-additive XC and NAKP con- 
tributions. Previous WFT-in-DFT implementations for 
open-shell systems have in practice neglected this differ- 
ence, effectively assuming that the spin polarization is 
localized within the WFT subsystem! 34 ! 53 ! In this study, 
we show that effects due to spin-dependent potentials 
are substantial and easily included via separate a and 
/3 embedding potentials. We develop approaches to uti- 
lize both restricted and unrestricted open-shell orbital 
formulations in WFT calculations. 



1. Open-shell DFT-in-DFT embedding 

We begin by considering an open-shell system for 
which the total electronic density is comprised of the a 
and P density of the two subsystems, pab = Pa + Pa + 
Pb + Pr- The effective potential for the unrestricted spin 
orbitals^ is 

v&\p%,pi,(%,&;r] = v^- a [pl,p{;r] +< mb (r)(ll) 

where v^' a [p\, p A ; r] is the standard KS effective poten- 
tial for the unrestricted (U)KS orbitals, and v° mh (r) is a 
spin-dependent embedding potential applied only to the 
a-spin electrons, 

Cb(r) = + Mpb; r] + <Mb= Pab; r] - 

vZ[pl,pi;r]+v^ d [pl,pt B ;r}. (12) 

The corresponding quantities for the /3-spin electrons are 
analogously defined. The total energy functional for the 



full open-shell system is then 

E[p AB ] = T s [p a A ] + T s [pg] + TT d [PA,Pg] + 
T s [ PA ]+T s [p^+Tr d [p A ,P B} + 
V ne [pAB] + J[pab] + E xc 

Ip1b,pU- (13) 

Separate OEP calculations are performed over the a and 
j3 spin-densities for the exact calculation of the NAKP, 
which allows for numerically exact unrestricted open- 
shell DFT embedding (U-DFT-in-DFT) P 

In practice, we solve for the unrestricted spin orbitals 
by adding the spin-dependent embedding potentials to 
the a and j3 Fock matrices. The a and /3 embedding 
potential can be written in the AO basis as 

v Lb = v nc + J[7B] + v« c [7ab, Tab] " 

vL^Til + vLibLTAB], ( 14 ) 

where £ <E {a, /?}, and the corresponding Fock matrices 
are 

f€,A in B _ j^A _|_ j[7a] + ^j^^ + ^ (15) 

The unrestricted spin orbitals for subsystem A are 
then obtained by separately diagonalizing f Q,A 111 B and 
f/3A 111 B in the usual way. 

Practical implementations for performing OEP calcu- 
lations using restricted open-shell orbitals have yet to 
be developed. We thus introduce a simple, approximate 
scheme for restricted open-shell DFT embedding (RO- 
DFT-in-DFT). In this approach, a U-DFT-in-DFT cal- 
culation is first performed, and the embedd ing poten- 
tials v™ mb and vf mb are constructed using Eq. |l4[ T hen, 
f/a,A 111 B and f^' A 111 B are constructed using Eq. |15[ and 
the usual RO approach is employed to obtain subsystem 
orbitals that are spatially identical for the a and j3 elec- 
trons. Specifically, f a,A 111 B is diagonalized to obtain a 
set of occupied a spin orbitals, {ft*}, and f^' A 111 B is 
then projected into the space of the occupied a spin or- 
bitals using 

F' A in B = c T f* A in B c, (16) 

where c is the matrix with columns comprised of the AO 
coefficients for {4>^ A }- Finally, the projected Fock ma- 
trix, f^' A 111 B , is diagonalized to obtain the set of RO 
orbitals, {</> A cc }, with the first N 13 ^ orbitals doubly oc- 
cupied and with orbitals A^' A + 1, . . . ,iV Q,A singly oc- 
cupied, where N a ' A and N@' A indicate the number of a 
and electrons in subsystem A. Although the second and 
third terms on the RHS of Eq. [15] are updated at each 
iteration of the RO-DFT-in-DFT calculation, we leave 
the embedding potentials unchanged to avoid perform- 
ing OEP calculations using restricted open-shell orbitals. 
The RO-DFT-in-DFT energy for the total density is cal- 
culated using Eq. [l3j 

Several different schemes have been proposed to calcu- 
late the embedding potential for ope n-shell subsystems 
while neglecting its spin- dependence! 34 ' 35 ' 53 ' These ap- 
proaches generally assume that interactions between the 
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subsystems can be described by a single embedding po- 
tential. For example, in systems with an even number 
of electrons, the embedding potential, v om b in Eq. [7j can 
be obtained assuming that each embedded subsystem is 
closed-shell, and then the open-shell subsystem is cal- 
culated using v£ mb = vf mb = v crab . In this approach, 
Eq, 
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is solved self-consistently while v" mb and v£ mb 
are held fixed, and the final DFT-in-DFT energy is calcu- 
lated using Eq. [13] This spin-independent description of 
the embedding potential can be used with either an unre- 
stricted treatment of the open-shell subsystem (U-DFT- 
in-DFT-CS) or a restricted treatment of the open-shell 
subsystem (RO-DFT-in-DFT-CS); we later employ the 
approach to compare with the previously described meth- 
ods (U-DFT-in-DFT and RO-DFT-in-DFT) that include 
spin-dependence in the embedding potential. 



2. Open-shell WFT-in-DFT embedding 

Unrestricted open-shell WFT-in-DFT (U-WFT-in- 
DFT) calculations are performed by first computing 
unrestricted Hartree-Fock (UHF) orbitals in the spin- 
dependent embedding potential, and then using these or- 
bitals for a post-HF WFT calculation. The a and f3 Fock 
matrices for the calculation of the UHF orbitals are 



= h A +J[ 7A ]+K«[ 7 i]+v 



cmb' 



(17) 



where £ € {a, /?}, K is the HF exchange matrix, and the 
embedding potentials (Eq. 14) are obtained from a U- 
DFT-in-DFT embedding calculation. The total energy 
is evaluated using 

P r^WFT DFT] _ pDFTt DFTi 

-frtot [Pa ) Pb J - ^ab [Pab I 



< FT br T ] 



E 
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«e{a,/3} 
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( r )PA DFT ( r ) dr 



"cmb 



r)pi WFT (r)dr|.(18) 



Restricted open-shell WFT-in-DFT (RO-WFT-in- 
DFT) calculations are performed by solving for restricted 
open-shell HF (ROHF) orbitals in the spin-dependent 
embedding potential. Just as in RO-DFT-in-DFT em- 
bedding, a U-DFT-in-DFT calculation is first performed, 
and the embedding potentials v" mb and v^ mb are con- 
structed using Eq. [14] Then, the Fock matrices in Eq. [IT] 
are constructed and the usual approach is employed to 
obtain RO orbitals; the second and third terms on the 
RHS of Eq. [T7j are updated at each iteration while the 
embedding potential is left unchanged. The ROHF or- 
bitals are used in the post-HF WFT calculation, and the 
total energy is then evaluated using Eq. [18] We note that 
for the RO- WFT-in-DFT energy calculation, the term 
e a FT [Pa FT }i is evaluated using the ROKS-DFT energy. 



B. Optimized effective potential 

As seen in Eqs.[5]and[6j DFT embedding requires com- 
putation of both the kinetic energy, T s [pab], and its func- 
tional derivative. However, since the explicit functional 
form for the kinetic energy is unknown, OEP methods 
are needed to obtain these terms exactly. 

The OEP is the local potential for which solution of 
the one-electron equations 



voEp(,rj 



(19) 



yields orbitals that correspond to a given target density 
while minimizing the non-interacting kinetic energy. A 
variety of methods for determining such potentials from 
an input target density have been developedP^HlH Cal- 
culations reported here employ the d irect optimization 
procedure developed by Wu and Yang! i n which the 
kinetic energy is obtained via the unconstrained maxi- 
mization 

T s [ PiD ]= max {W s [y dct ,v OE p(r)]}, (20) 

»OEp(r) 

where 



W / s [*dct,«OEp(r)] = 2 



(poep(i\) - Pin(r)) uoEp(r)dr 
Cl|V«A(r)|| 2 , (21) 



and 



^OEp(r) = v c £[p hl ;r] + u A (r). 



(22) 



In these equations, v\(r) — X)t^t3*( r ); {9t( r )} comprise 
an auxiliary basis set for the potential, b t are the cor- 
responding expansion coefficients, and £ is a regulariza- 
tion parameter.^ Maximization of W s utilizes the New- 
ton method for optimization with a back-tracking line 
search in the expansion coefficients,^ 



b «+i) 



(23) 



where i is the iteration number, H and g are the Hessian 
and gradient of W s , respectively, and r € [0,1] is the 
step-size in the line search. 

In practice, to obtain the embedding potential, we do 
not explicitly calculate the NAKP for each subsystem. 
Instead, for closed shell subsystems, we directly update 
the embedding potential (Eq. [4]) at each iteration of the 
KSCED equations using^ 



^emb 



(r) 



, (») 
cmb 



(r) - 0v x {i 



(24) 



where 9 G [0,1] is a damping coefficient. By construction, 
the embedding potential for each subsystem is identical 
at every iteration. The KSCED equations are initialized 
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using u^ b (r) = 0, such that the initial guess for the 
NAKP for subsystem A exactly cancels the remaining 
terms in Eq. [4] (and the initial guess for the NAKP for 
subsystem B likewise cancels the corresponding terms). 
Upon convergence of the KSCED equations, v\(r) = 0, 
[pab I r] = vf£[p KS ;r], and v emh (r) = ^ b (r). En- 



forcing the embedding potential to be identical for all 
subsystems l eads t o a unique partitioning of the subsys- 
tem densitiesPfiEU 

For open-shell calculations, we similarly update the 
spin-dependent embedding potential (Eq. 12); the OEP 
obtained for a given spin density is 



^OEp( r ) = 



KS,{r/ c 

'eg m 



+ PB 



•§),(/£ + pg);r]+«l(r), (25) 



and as in Eq. |24| v\(r) is used to update wf mb (r) a * eacn 
iteration. 

Finally, we note that XC functionals that include a 
fraction of the exact exchange can be employed in DFT 
embedding via the OEP calculation. The HF exchange 
matrix, K, is evaluated using 7oePj the OEP density 
matrix in the AO basis. For DFT-in-DFT embedding, 
the exchange energy is calculated using 

Ex [70EP,7in] = -tr(7oEpK[7oEp]) 

+tr ((7oep - 7m)K[ 70 Ep]) , (26) 

where the second term on the RHS corrects the exchange 
energy for small numerical differences between 7oep and 
7i n . For calculations on the low-spin state of the hex- 
aaquairon(II) cation, this correction is found to be as 
large as 20 kcal/mol; however, the correction is not re- 
quired for the evaluation of the WFT-in-DFT energy 



(Eq. 10), since ^ab T [Pab T ] ^ s obtained directly from a 



KS-DFT calculation on the full system. 



C. Orbital-occupation freezing 

For Ws to be a concave function of w OEP (r), it is 
necessary^ that the orbitals used to construct poep in 
Eq. [21] correspond to the lowest eigenvalues of Eq. [T9| 
However, this can be problematic for systems with small 
energy differences between the occupied and virtual or- 
bitals, where small changes in uoEp(r) can alter the rel- 
ative ordering of the orbitals. 

To illustrate this issue, Fig. [T] shows the line search 
for an illustrative Newton step in an OEP calculation 
for the low-spin hexaaquairon(II) cation. W s is plotted 
as a function of r, where r is the step-size in Eq. |23| 
For any step-size larger than r = 0.38 in this case, the 
orbital occupancy changes from one in which only t2 g -like 
d orbitals are occupied to one in which e g -like d orbitals 
are occupied. In traditional back-tracking line searches, 
any step which increases W s would be accepted, including 
the t = 0.5 step indicated with the red arrow. However, 
this step is problematic since the Hessian and gradient of 
W s for the next Newton step would be evaluated using a 
density that corresponds to the wrong orbitals. The net 
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FIG. 1. An illustrative Newton step in the OEP calcu- 
lation for the low-spin hexaaquairon(II) cation, performed 
with (black) and without (red) the orbital-occupation-freezing 
technique. The technique ensures that correct orbitals remain 
occupied throughout the maximization of W s . See text for de- 
tails. 



results are poor convergence and incorrect solutions for 
the OEP. 

We introduce a simple method to alleviate this prob- 
lem by modifying the back-tracking line search. Ref- 
erence (r=0) orbitals are computed from Eq. 19 us- 
ing woEp(r) = w^ S [pab; r], and for any proposed step- 
size t, the corresponding orbitals are computed using 



^OEp(r) 



KS 
"eff 



[pab! 1- ] + v \( r )- The proposed step is 



rejected if the overlap between these two sets of orbitals 
is less than 0.5, regardless of the change in W s ; otherwise, 
it is subjected to the usual criteria of the back-tracking 
line search algorithm. Upon rejection, the step-size r is 
reduced by a factor of 2. This technique ensures that the 
correct orbitals remain occupied throughout the maxi- 
mization of W s - In Fig. [l] the proposed step indicated 
by the red arrow is rejected, whereas the proposed shorter 
step indicated by the black arrow is accepted; not only is 
the value of Ws increased, but the correct orbitals remain 
occupied. By utilizing this technique, we found that the 
maximization of W s typically requires less than 20 New- 
ton steps for the low spin state of the hexaaquairon(II) 
cation, whereas the optimization failed to converge with- 
out the use of orbital-occupation freezing. 



D. Computational Details 

The DFT embedding methods employed here are all 
implemented in the development version of the Molpro 
software packaged All calculations employ the super- 
molecular basis set convention, in which the molecular 
orbitals for each subsystem are described in the AO ba- 
sis for the full system.^ Calculations on the ethylene- 
propylene dimer use the aug-cc-pVTZ orbital basis set 
for the carbon atoms and the aug-cc-pVDZ orbital basis 
set for the hydrogen atoms. Calculations on the hex- 
aaquairon(II) cation use the aug-cc-pVTZ orbital basis 
set for the iron atom and the aug-cc-pVDZ orbital basis 
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set for the hydrogen and oxygen atoms. For the aux- 
iliary basis set used in the OEP calculations, we em- 
ploy atom-centered Gaussian basis functions (<?j(r) = 
N t e~ x±r , where N t is the normalization constant) for 
which the coefficient A t assumes values of 2™, where 
n = n min , n min + 2,..., n max - 2, n max . Calculations on 
the ethylene-propylene dimer employ the basis set for 
which the s-type functions for the carbon and hydrogen 
atoms span {n m i n ,n max } = {-4, 4}, and the p-type func- 
tions for the carbon and hydrogen atoms span {-2, 2}. 
Calculations for the hexaaquairon(II) cation employ the 
basis set for which the s-type functions for the iron atom 
span {-4, 6}, the p-type functions for the iron atom span 
{-4, 6}, the d-type functions for the iron atom span {-2, 
2}, the s-type functions for the oxygen atoms span {-4, 
6}, the p-type functions for the oxygen atoms span {-2, 
4}, the s-type functions for the hydrogen atoms span {- 
4, 4}, and the p-type functions for the hydrogen atoms 
span {-2, 2}. For all systems, the finite auxiliary ba- 
sis set for the OEP calculations was confirmed to intro- 
duce a difference of less than 1 kcal/mol between the 
total energy computed using KS-DFT and either closed- 
shell or unrestricted open-shell DFT-in-DFT embedding. 
The regularization parameter used in the OEP calcula- 
tions is set to C = 10~ 3 ; smaller values were tested on 
the ethylene-propylene dimer and the hexaaquairon(II) 
cation and were found to have only a small (0(/iHartree)) 
effect on the total DFT-in-DFT energy. 

The KSCED equations are initialized with subsystem 
densities comprised of the superposition of HF atomic 
densities and with u em b(r) = 0; different initial guesses 
for the embedding potential were tested on the hex- 
aaquairon(II) cation and were found to yield similar final 
embedding potentials with only small (O(10 ^Hartree)) 
changes in the total DFT-in-DFT energy. 
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FIG. 2. WFT-in-DFT embedding for the ethylene-propylene 
dimer. (a) The ethylene-propylene dissociation curve, ob- 
tained using CCSD(T)-in-B3LYP (red) and KS-DFT with 
PBE (green), B3LYP (orange), B-LYP (blue) and B88-P86 
(cyan) for the XC functional. Also included are the reference 
CCSD(T) results (black), which are graphically indistinguish- 
able from the CCSD(T)-in-B3LYP results. The curves are 
vertically shifted to align at infinite separation, (b) Isosur- 
face plots indicate the subsystem partitioning for the ethene- 
propene dimer calculations. The red isosurface indicates the 
density of the 32 electrons associated with the C2H4-C2H3- 
moiety, and the blue isosurface indicates the density of the 
8 electrons associated with the -CH3 moiety. The isosurface 
plot corresponds to an electronic density of 0.05 a.u. 



IV. RESULTS 

A. The Ethylene-Propylene Dimer: WFT-in-DFT 
Embedding 

The ethylene-propylene dimer is a prototypical sys- 
tem for which quantum embedding methods, such as 
QM/MM or ONIOM, may be employed, ft exhibits a 
weak 7T — 7T interaction that is difficult to address with 
conventional KS-DFT methods, while also exhibiting a 
spectator -CH 3 moiety that contributes little to the in- 
teraction energy while substantially increasing the cost of 
the high-level calculation. However, unlike the QM/MM 
treatment of subsystems, the interactions between the 
7r — 7r system and the -CH3 moiety can be treated seam- 
lessly using WFT-in-DFT embedding, as is now demon- 
strated. 

Fig. [2ja) presents the ethylene-propylene dimer dis- 
sociation curve plotted as a function of the distance be- 
tween the ethylene and propylene ir bonds, with the equi- 
librium dimer geometry obtained via minimization at the 



MP2 level of theory. Other geometries along the curve 
are obtained by displacing the two molecules along the 
vector formed between the midpoints of the two C=C 
bonds, while fixing all other internal coordinates. The 
relative energies are plotted by aligning each curve at in- 
finite separation. The full CCSD(T) calculation (black) 
shows a binding energy of 2.0 kcal/mol. KS-DFT cal- 
culati ons us ing the PBlPMl (green), B3LYP^ (orange), 
B-LYPSMH (blue) and B88-P8#Ha ( cy an) XC function- 
als illustrate the difficulty in describing dispersion inter- 
actions using KS-DFT. The PBE functional underesti- 
mates the binding energy by 1.3 kcal/mol, while the rest 
of the XC functionals fail to capture any of the attractive 
interactions. 

Finally, the red curve in Fig. [2]ja) presents the results 
of WFT-in-DFT embedding, using a subsystem parti- 
tioning in which the 32 electrons associated with the 
7r system (C2H4-C2H3-, red in Fig. [2jb)) are treated 
at the WFT level of theory and the remaining 8 elec- 
trons in the -CH3 moiety are treated at the DFT level 
of theory. We employ CCSD(T) for the WFT and the 
B3LYP XC functional for the DFT (i.e., CCSD(T)-in- 



r 



B3LYP). Fig.[2](a) shows excellent agreement between the 
CCSD(T) (black) and CCSD(T)-in-B3LYP (red) calcu- 
lations; these curves, which are graphically indistinguish- 
able, differ by less than 0.10 kcal/mol through the entire 
range of distances. We have confirmed that this level 
of accuracy is maintained with different XC functionals 
used for the DFT; specifically, CCSD(T)-in-(B-LYP) en- 
ergies differ from the CCSD(T) results by less than 0.20 
kcal/mol throughout the entire curve. These results il- 
lustrate that WFT-in-DFT embedding can be used to 
systematically improve DFT results and to avoid embed- 
ding errors while partitioning across covalent bonds. 



B. The hexaaquairon(ll) cation 

We now present DFT-in-DFT and WFT-in-DFT cal- 
culations for the high-spin [ 5 T 2g : (t 2g ) 4 (e g ) 2 ] and low- 
spin [ 1 A lg : (t 2 g) 6 (e g ) ] states of the hexaaquairon(Il) 
cation, a system that presents challenges due to the pres- 
ence of low-lying unoccupied orbitals, the important role 
of unpaired electrons, and the relatively large number of 
electrons (84 e~) in the full system. First, we test the ac- 
curacy of DFT-in-DFT embedding for the various treat- 
ments of the open-shell embedding potential described 
earlier. We then employ WFT-in-DFT calculations to 
investigate the low-spin/high-spin energy splitting and 
the ligation energy for this transition metal complex. 



1. DFT-in-DFT embedding 

Fig. |3ja) presents the potential energy curve for the 
simultaneous dissociation of all six H 2 ligands of the 
hexaaquairon(II) cation, plotted as a function of the 
average iron-oxygen distance. The equilibrium geome- 
tries for the low-spin [ 1 A ig : (t 2g ) 6 (e g ) ] and high-spin 
[ 5 T 2g : (t 2g ) 4 (e g ) 2 ] states are obtained using KS-DFT 
energy minimization with the B3LYP XC functional; all 
other geometries are obtained by uniformly stretching the 
iron-oxygen distances in the complex, keeping all other 
internal coordinates unchanged. All KS-DFT and DFT- 
in-DFT embedding results reported in this section are 
obtained using the B3LYP XC functional. The curves 
in the main panel of Fig. [3]ja) are vertically shifted to 
share a common minimum value; they are not horizon- 
tally shifted. The high-spin state is lower in energy and 
exhibits a longer average iron-oxygen distance than the 
low-spin state. 

We perform DFT-in-DFT embedding using a subsys- 
tem partitioning in which the 24 electrons associated 
with the iron center comprise one subsystem (red in 
Fig.[3^b)) and the remaining 60 electrons associated with 
the six water ligands comprise a second subsystem (blue 
in Fig. [3^b)). For the low-spin state, Fig. [3](a) demon- 
strates good numerical agreement between DFT-in-DFT 
(red) and KS-DFT (black) ; the relative energies differ by 
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FIG. 3. DFT-in-DFT embedding for the hexaaquairon(II) 
cation, (a) The potential energy curve for the simultane- 
ous dissociation of the six H2O ligands. All curves in the 
main panel are vertically shifted to share a common minimum 
energy; they are not horizontally shifted. The dissociation 
curves for the low-spin ( 1 A ig ) state obtained using KS-DFT 
(black) and DFT-in-DFT (red) are graphically indistinguish- 
able. The dissociation curves or the high-spin ( T2 g ) state 
obtained using UKS-DFT (blue), U-DFT-in-DFT (green), 
ROKS-DFT (magenta), and RO-DFT-in-DFT (orange) are 
likewise graphically indistinguishable. The inset shows these 
four high-spin potential energy curves, with each curve ver- 
tically shifted only by the UKS-DFT minimum energy of 
— 1721.693423 Hartree. The dashed black dissociation curve 
in the main panel is obtained using the RO-DFT-in-DFT-CS 
method, which neglects spin-dependence in the embedding 
potential, (b) Isosurface plots indicate the subsystem par- 
titioning for the hexaaquairon(II) cation. The red isosurface 
indicates the density of the 24 electrons associated with the Fe 
atom, and the blue isosurface indicates the density of the 60 
electrons associated with the six H2O ligands. The isosurface 
plot corresponds to an electronic density of 0.05 a.u. 



less than 0.6 kcal/mol throughout the range of reported 
internuclear distances. 

For the high-spin state of the hexaaquairon(II) cation, 
Fig. ga) shows that the UKS-DFT and ROKS-DFT 
methods are in good agreement with each other, as well 
as with the corresponding U-DFT-in-DFT and RO-DFT- 
in-DFT embedding approaches described in Sec. Ill A 
1. The U-DFT-in-DFT calculation accurately reproduces 
the relative energies obtained from UKS-DFT to within 
0.4 kcal/mol throughout the attractive branch of the 
curve and to within 0.8 kcal/mol at shorter distances. 
The RO-DFT-in-DFT calculation reproduces the relative 



energy obtained from ROKS-DFT to within 1.0 kcal/mol 
throughout the attractive branch of the curve and to 
within 2.2 kcal/mol at shorter distances. 

The inset of Fig. [3^a) shows the various potential en- 
ergy curves computed for the high-spin state of the hex- 
aaquairon(II) cation, with each curve vertically shifted by 
only the UKS-DFT minimum energy. This inset demon- 
strates relatively small differences in the total energies 
computed with the various embedding and open-shell 
treatments. 

Finally, the dashed black curve in Fig. [3]Ja) demon- 
strates the importance of including spin-dependence in 
the embedding potential. This curve corresponds to the 
RO-DFT-in-DFT-CS treatment of the high-spin state of 
the hexaaquairon(II) cation described in Sec. Ill A 1. It 
exhibits large relative errors (over 70 kcal/mol) compared 
to the other treatments of the high-spin state of the hex- 
aaquairon(II) cation, as well as qualitatively incorrectly 
shortening of the equilibrium internuclear distance. Al- 
though this approximation is expected to be more reli- 
able for systems in which the spin-density is strongly lo- 
calized with a single subsystem, the result demonstrates 
that substantial errors can emerge due to the neglect of 
spin-dependence in the embedding potential. 



2. WFT-in-DFT Embedding 

We now consider WFT-in-DFT embedding for the 
hexaaquairon(II) cation, employing the same subsystem 
partitioning as in the DFT-in-DFT embedding calcu- 
lations (Fig. [3^b)). The hexaaquairon(II) cation is a 
benchmark system for spin splittings in transition metal 
complexesPS We initially discuss results for MP2 embed- 
ding to compare the U- WFT-in-DFT and RO- WFT-in- 
DFT approaches, and we then present results obtained 
using CCSD(T) embedding. 

Fig. [2] presents results for the low-spin/high-spin en- 
ergy difference (AE hH ) obtained using MP2, KS-DFT, 
and MP2-in-DFT embedding; detailed values are re- 
ported in Table [TJ For KS-DFT calculations of A^lh, 
the energy for the high-spin state of the hexaaquairon(II) 
cation was obtained at the UKS-DFT level of theory. The 
WFT-in-DFT embedding energy for the low-spin state 
of the hexaaquairon(II) cation is obtained using closed- 
shell WFT-in-DFT (Sec. II B), while the high-spin state 
is treated using cither U- WFT-in-DFT or RO- WFT-in- 
DFT (Sec. Ill A 2). The KS-DFT results (red in Fig. Q 
exhibit strong dependence on the XC functional, with 
hybrid functionals underestimating AE'lh to a somewhat 
lesser degree than the semi- local functionals. 

Fig. H| clearly illustrates that the RO-MP2-in-DFT re- 
sults (blue) are in better agreement with the full MP2 
calculation than the corresponding U-MP2-in-DFT re- 
sults (green), particularly for semi- local XC functionals. 
Removal of spin-contamination in the WFT calculation 
reduces the energy of the high-spin state RO- WFT-in- 
DFT calculation with respect to that obtained using U- 



TABLE I. High-spin/low-spin splitting energies in cm 1 for 
the hexaaquairon(II) cation obtained using KS-DFT, U-MP2- 
in-DFT, and RO-MP2-in-DFT with a range of different XC 
functionals. ft b 



Functional 


KS-DFT 


U-MP2-in-DFT 


RO-MP2-in-DFT 


B-LYP 


7828 


12604 


15294 


PBE 


9479 


11079 


13395 


PW91 


8593 


10924 


13201 


B3LYP 


11206 


14387 


14703 


PBEO 


14154 


13812 


13979 



a RO-MP2 yields 16439 cm" 1 . 
b U-MP2 yields 17396 cm" 1 . 



WFT-in-DFT. 

Another important observation from Fig. [3] is that 
the dependence of Ai?LH on the DFT XC functional 
is greatly reduced in the embedding calculation, even 
though only the single transition metal atom is treated 
at the WFT level. The spread of values obtained at the 
KS-DFT level of theory is over 6000 cm -1 , which is re- 
duced by a factor of 3 in the RO-MP2-in-DFT embedding 
calculations. 
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FIG. 4. MP2-in-DFT embedding for the hexaaquairon(II) 
cation. High-spin/low-spin splitting energies obtained using 
KS-DFT (red,circles), U-MP2-in-DFT (green, squares), and 
RO-MP2-in-DFT (blue,triangles) with a ran ge of different 
XC functionals that include B-LYPF^ PBEp^Ml PW91, 1111 
B3LYP|^and PBEO The black line indicates the reference 
value of 16439 cm - obtained at the RO-MP2 level of theory; 
U-MP2 yields a value of 17396 cm -1 . 



Fig. |5ja) presents calculations of the low-spin/high- 
spin splitting obtained using WFT-in-DFT calculations 
at the RO-CCSD(T)-in-DFT level of theory; detailed val- 
ues are reported in Table [Til F° r the reference calculation 
obtained at the full RO-CCSD(T) level of theory,^ no T2 
amplitudes were found to exceed 0.05, indicating that a 
single-reference description of the wavefunction is ade- 
quate. The general trend for the RO-CCSD(T)-in-DFT 
calculations is consistent with the results obtained from 
RO-MP2-in-DFT. It is again seen that the dependence of 
Ai/LH on the XC functional is substantially reduced us- 
ing RO-CCSD(T)-in-DFT embedding, and the accuracy 
of the KS-DFT results are generally improved by treat- 
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ing the transition metal atom at the WFT level. For 
this system, the embedded RO-CCSD(T) calculation in- 
volves correlating significantly fewer electrons than the 
full RO-CCSD(T) calculation, and we found that the 
WFT step in the RO-CCSD(T)-in-DFT calculation re- 
quired approximately 50 times less wall-clock time than 
the full RO-CCSD(T) calculation. 

Fig. [5jb) shows that the LDA functiona l 68 * 69 1 presents 
an interesting outlier compared to the other results in 
Fig.ga). Unlike the semi-local and hybrid functionals, 
RO-CCSD(T)-in-LDA calculations do not exhibit a sig- 
nificant improvement with respect to the corresponding 
KS-DFT result. We now show that this anomalous result 
arises from a density-based error in the LDA functional. 

Fig. [5jc) and Fig. [5fd) present the charge on the Fe 
atom from a Mulliken population analysis for the low- 
spin state of the hexaaquairon(II) cation. Fig.^c) shows 
that the semi-local and hybrid functionals all yield a sim- 
ilar charge for the Fe atom, which is very close to that of 
the full (relaxed) CCSD density. In contrast, Fig. [5jcL) 
reveals the LDA functional significantly underestimates 
the Fe atomic charge, which indicates a significant error 
in the calculation of the ground state density. Although 
the use of embedded WFT can be expected to overcome 
the error in the contribution to the spin-splitting energy 
due to the LDA functional, it can not overcome this error 
in the actual ground state density due to LDA. 

To confirm this interpretation, we show that removing 
the error in the LDA density leads to improved WFT-in- 
DFT estimates for the spin-splitting energy, even if the 
LDA functional is still employed for the DFT contribu- 
tions to the energy. In Fig.[5jb), the B-LYP+LDA result 
for WFT-in-DFT embedding (blue, triangle) is obtained 
by (i) calculating the embedding potential and the sub- 
system densities using the B-LYP XC functional, (ii) per- 
forming the embedded WFT calculation at the CCSD(T) 
level, and (Hi) using the LDA functional and CCSD(T) 
to evaluate the respective DFT and WFT contributions 
to the total energy in Eq. [l8j The corresponding B- 
LYP+LDA result for KS-DFT (red, circle) is obtained 
by calculating the total density using KS-DFT with the 
B-LYP XC functional and then using the LDA functional 
to evaluate the KS-DFT energy. As is seen in Fig. [5jd), 
the B-LYP treatment of the subsystem densities leads to 
the expected partial charge for the Fe atom; it avoids 
the error in the electronic density that is introduced us- 
ing LDA. However, the spin-splitting energy obtained us- 
ing the B-LYP+LDA result for KS-DFT is essentially no 
better than that obtained using KS-DFT with the LDA 
functional (Fig[5jb)), indicating that simply correcting 
the LDA error in the density is not enough to avoid 
the LDA error in the energies. Finally, Fig. [5jR) shows 
that the B-LYP+LDA result for WFT-in-DFT does ex- 
hibit a substantial improvement over the corresponding 
KS-DFT result; this confirms that WFT embedding is 
able to overcome energy-based errors due to the DFT 
XC functional, although it is less effective at overcoming 
density-based errors due to the DFT XC functional. 
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FIG. 5. CCSD(T)-in-DFT embedding for the hexaaqua- 
iron(II) cation. (a,b) High-spin/low-spin splitting energies 
obtained using KS-DFT (red,circles) and RO-CCSD(T)-in- 
DFT (blue,triangles) with a range of different XC function- 
als. The B-LYP+LDA result is obtained using the B-LYP 
XC functional for the density calculation and the LDA XC 
functional for the energy calculation, as is described in the 
text. The black line indicates the reference value of 14149 
cm" 1 obtained at the RO-CCSD(T) level of theory. (c,d) The 
charge on the Fe atom is obtained using the Mulliken popula- 
tion analysis of the KS-DFT calculation with each functional. 
The relaxed CCSD density, indicated by the black line, has 
an Fe atomic charge of 2.56. 



TABLE II. High-spin/low-spin splitting energies in cm 1 for 
the hexaaquairon(II) cation obtained using KS-DFT and RO- 
CCSD(T)-in-DFT with a range of different XC functionals. Q 



Functional 


KS-DFT 


RO-CCSD(T)-in-DFT 


B-LYP 


7828 


12554 


PBE 


9479 


11238 


PW91 


8593 


10712 


B3LYP 


11206 


12634 


PBEO 


14154 


12912 



RO-CCSD(T) yields 14149 cm" 



Although we have shown that WFT-in-DFT embed- 
ding with the subsystem partitioning shown in Fig. (3^b) 
generally leads to improved estimates for the low- 
spin/high-spin splitting energy over KS-DFT, the same 
does not hold true for calculated ligation energies of the 
hexaaquairon(II) cation. Ligation energies calculated us- 
ing RO-CCSD(T)-in-DFT embedding are essentially un- 
changed from those obtained using KS-DFT with the 
corresponding XC functional; indeed, the mean abso- 
lute difference between the computed WFT-in-DFT and 
KS-DFT ligation energy is only 0.6 kcal/mol per lig- 
and across the set of functionals that includes LDA, B- 
LYP, PBE, PW91, B3LYP, and PBEO. Unlike the spin- 
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splitting energy, which is highly sensitive to the elec- 
tronic structure of the Fe atom and is thus impacted 
by the WFT subsystem description, the ligation energy 
is dominated by interactions between the Fe atom and 
the water ligands; these inter-subsystem interactions are 
still treated essentially at the DFT level in WFT-in-DFT 
embedding. An improved description for the ligation en- 
ergy could be obtained by simply expanding the number 
of electrons that are treated at the WFT level of theory, 
or by including two-body correlation corrections through 
an embedded many-body expansion description of the 
systemP3 



CONCLUSIONS 

In this work, we have introduced and demonstrated im- 
proved methods for the implementation of WFT-in-DFT 
calculations for open-shell systems and systems with 
low-lying virtual orbitals. A simple orbital-occupation- 
freezing technique is introduced to enable robust OEP 
calculations on systems with small HOMO-LUMO gaps, 
leading to accurate DFT-in-DFT and WFT-in-DFT em- 
bedding calculations on transition-metal complexes. Fur- 
themore, the use of spin-dependent embedding potentials 
is shown to preserve the accuracy of open-shell DFT-in- 
DFT calculations in both the restricted and unrestricted 
orbital formulations, whereas neglect of the spin polariza- 
tion leads to significant errors in both computed energies 
and geometries. WFT-in-DFT calculations on the hex- 
aaquairon(II) cation reveal that the treatment of only the 
single transition metal atom leads to significant improve- 
ments in the accuracy of calculated spin-splittings, as 
well as marked reduction in the dependence of results on 
the DFT XC functional. Taken together, the exact em- 
bedding techniques reported and demonstrated here offer 
a promising approach to the robust treatment of systems 
for which the accuracy of WFT is required but for which 
the cost of the full WFT calculation is not feasible. 



ACKNOWLEDGEMENTS 

This work is supported in part by the Air Force Of- 
fice of Scientific Research (FA9550-11-1-0288) and the 
U. S. Army Research Laboratory and the U. S. Army 
Research Office (W911NF-10-1-0202). TFM and FRM 
also gratefully acknowledge network funding from the 
NSF (CHE-1057112) and EPSRC (EP/J012742/1), re- 
spectively. Computational resources were provided by 
the National Energy Research Scientific Computing Cen- 
ter, which is supported by the Office of Science of the US 
Department of Energy under Contract No. DE-AC02- 
05CH11231. 

1 A. Warshel and M. Levitt, J. Mol. Biol. 103, 227 (1976). 

2 P. Sherwood, A. H. de Vries, S. J. Collins, S. P. Greatbanks, 

N. A. Burton, M. A. Vincent, and I. H. Hillier, Farad. Discuss. 

106, 79 (1997). 



3 J. L. Gao, P. Amara, C. Alhambra, and M. J. Field, J. Phys. 

Chem. A. 102, 4714 (1998). 
4 H. Lin and D. G. Truhlar, Theor. Chem. Ace. 117, 185 (2007). 
5 H. M. Senn and W. Thiel, Angew. Chem. Int. Ed. 48, 1198 

(2009). 

6 L. Hu, P. Soderhjelm, and U. Ryde, J. Chem. Theory. Comput. 
7, 761 (2011). 

7 S. Dapprich, I. Komaromi, K. S. Byun, K. Morokuma, and 

M. J. Frisch, J. Mol. Struct. THEOCHEM 461-462, 1 (1999). 
8 F. Maseras and K. Morokuma, J. Comp. Chem. 16, 1170 (1995). 
9 K. Kitaura, E. Ikco, T. Asada, T. Nakano, and M. Uebayasi, 

Chem. Phys. Lett. 313, 701 (1999). 
10 D. G. Fedorov and K. Kitaura, J. Chem. Phys. 120, 6832 (2004). 
n D. G. Fedorov and K. Kitaura, J. Phys. Chem. A. Ill, 6904 

(2007) . 

12 M. J. Field, P. A. Bash, and M. Karplus, J. Comput. Chem. 11, 
700 (1990). 

13 U. C. Singh and P. A. Kollman, J. Comput. Chem. 7, 718 (1986). 
14 Y. Zhang, T. S. Lee, and W. Yang, J. Chem. Phys. 110, 46 
(1999). 

15 Y. Zhang, J. Chem. Phys. 122, 024114 (2005). 

16 G. Chalasinski, M. M. Szczesniak, P. Cieplak, and S. Schcincr, 

J. Chem. Phys. 94, 2873 (1991). 
17 E. B. Kadossov, K. J. Gaskell, and M. A. Langell, J. Comput. 

Chem. 28, 1240 (2007). 
18 G. Senatore and K. R. Subbaswamy, Phys. Rev. B. 34, 5754 

(1986). 

19 M. D. Johnson, K. R. Subbaswamy and G. Senatore, Phys. Rev. 

B. 36, 9202 (1987). 
20 P. Cortona, Phys. Rev. B 44, 8454 (1991). 

21 T. A. Wesolowski and A. Warshel, J. Phys. Chem. 97, 8050 
(1993). 

22 T. A. Wesolowski, in Computational Chemistry: Reviews of Cur- 
rent Trends - Vol. 10, pp. 1-82 (World Scientific, Singapore, 
2006). 

23 N. Govind, Y. A. Wang, A. J. R. da Silva, and E. A. Carter, 

Chem. Phys. Lett. 295, 129 (1998). 
24 P. Huang and E. A. Carter, J. Chem. Phys. 135, 194104 (2011). 
25 P. Elliott, K. Burke, M. H. Cohen, and A. Wasserman, Phys. 

Rev. A 82, 024501 (2010). 
26 J. Nafziger, Q. Wu, and A. Wasserman, J. Chem. Phys. 135, 

234101 (2011). 

27 J. D. Goodpaster, N. Ananth, F. R. Manby, and T. F. Miller, 
III, J. Chem. Phys. 133, 084103 (2010). 

28 J. D. Goodpaster, T. A. Barnes, and T. F. Miller, III, J. Chem. 
Phys. 134, 164108 (2011). 

29 F. R. Manby, M. Stella, J. D. Goodpaster, and T. F. Miller, III, 
J. Chem. Theory Comput. 8, 2564 (2012). 

30 O. Roncero, M. P. de Lara-Castells, P. Villarreal, F. Flores, J. Or- 
tega, M. Paniagua, and A. Aguado, J. Chem. Phys. 129, 184104 

(2008) . 

31 0. Roncero, A. Zanchet, P. Villarreal, and A. Aguado, J. Chem. 

Phys. 131, 234110 (2009). 
32 S. Fux, K. Kiewish, C. R. Jacob, J. Neugebauer, and M. Rciher, 

Chem. Phys. Lett. 461, 353 (2008). 
33 S. Fux, C. R. Jacob, J. Neugebauer, L. Visscher, and M. Reiher, 

J. Chem. Phys. 132, 164101 (2010). 
34 A. S. P. Gomes, C. R. Jacob, and L. Visscher, Phys. Chem. 

Chem. Phys. 10, 5353 (2008). 
35 C. Huang, M. Pavone, and E. A. Carter, J. Chem. Phys. 134, 

154110 (2011). 

36 A. Solovyeva, M. Pavanello, and J Neugebauer, J. Chem. Phys. 

136, 194104 (2012). 
37 C. R. Jacob, J. Neugebauer, and L. Visscher, J. Comput. Chem. 

29, 1011 (2008). 

38 M. Iannuzzi, B. Kirchner, and J. Hutter, Chem. Phys. Lett. 
421, 16 (2006). 

39 J. W. Kaminski, S. Gusarov, T. A. Wesolowski, and A. Ko- 
valenko, J. Phys. Chem. A. 114, 6082 (2010). 



11 



40 L. Rajchcl, P. S. Zuchowski, M. M. Szczgsniak, and 
G. Chalasinski, Chem. Phys. Lett. 486, 160 (2010). 

41 A. W. Gotz, S. M. Beyhan, and L. Visscher, J. Chem. Theory 
Comput. 5, 3161 (2009). 

42 A. S. P. Gomes and C. R. Jacob, Armu. Rep. Prog. Chem., Sect. 
C: Phys. Chem. 108, 222 (2012). 

43 Q. Wu and W. Yang, J. Phys. Chem. 118, 2498 (2003). 

44 F. A. Bulat, T. Heaton-Burgess, A. J. Cohen, and W. Yang, J. 

Chem. Phys. 127, 174101 (2007). 
45 Q. Wu, P. W. Ayers, and Y. Zhang, J. Chem. Phys. 131, 164112 

(2009). 

46 C. R. Jacob, J. Chem. Phys. 135, 244102 (2011). 
47 Q. S. Zhao and R. G Parr, Phys. Rev. A 46, 2337 (1992). 
4S Q. S. Zhao and R. G Parr, J. Chem. Phys. 98, 543 (1993). 
49 Q. S. Zhao, R. C. Morrison, and R. G Parr, Phys. Rev. A 50, 
2138 (1994). 

50 S. Sherifzadeh, P. Huang, and E. A. Carter, Chem. Phys. Lett. 
470, 347 (2009). 

51 T. Kliiner, N. Govind, Y. A. Yang, and E. A. Carter, J. Chem. 

Phys. 116, 42 (2002). 
52 T. A. Wesolowski, Phys. Rev. A 77, 012504 (2008). 
53 P. Huang and E. A. Carter, J. Chem. Phys. 125, 084102 (2006). 
54 S. Sherifzadeh, P. Huang, and E. A. Carter, J. Phys.: Condcns. 

Matter 21, 355501 (2009). 



55 A. Mordecai, in Nonlinear Programming: Analysis and Methods, 

(Dover Publishing, 2003). 
56 H.-J. Werner, P. J. Knowles, R. Lindh, F. R. Manby, M. Schiitz 

et al. Molpro, version 2012.1; Cardiff University: Cardiff, 

U. K.; Universitt Stuttgart: Stuttgart, Germany, 2012. See 

www.molpro.net. 
57 T. A. Wesolowski, J. Chem. Phys. 106, 8516 (1997). 
58 J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 

3865 (1996). 

59 J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 78, 
1396 (1997). 

60 A. D. Becke, J. Chem. Phys. 98, 5648 (1993). 
61 A. D. Becke, Phys. Rev. A 38, 3098 (1988). 

62 C. Lee, W. Yang, and R. G. Parr, Phys. Rev. B 37, 785 (1988). 

63 J. P. Perdew, Phys. Rev. B 33, 8822 (1986). 

64 J. P. Perdew, J. A. Chcvary, S. H. Vosko, K. A Jackson, 

M. R. Pederson, D. J. Singh, and C. Fiolhais, Phys. Rev. B. 

46, 6671 (1992). 
65 C. Adamo and V. Barone, J. Chem. Phys. 110, 6158 (1999). 
66 A. Fouqueau, S. Mer, M. E. Casida, L. M. L. Daku, A. Hauser, 

T. Mincva, and F. Neese, J. Chem. Phys. 120, 9473 (2004). 
67 P. J. Knowles, C. Hampel, and H.-J. Werner, J. Chem. Phys. 

99, 5219 (1993). 
68 J. C. Slater, Phys. Rev. 81, 385-390 (1951). 

69 S. H. Vosko, L. Wilk, and M. Nusair, Can. J. Phys. 58, 1200- 
1211 (1980). 



